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Abstract 

When applied to a single nucleon, nuclear energy density functionals may yield a non-vanishing 
internal energy thus implying that the nucleon is interacting with itself. It is shown how to avoid 
this unphysical feature for semi-local phenomenological functionals containing all possible bilinear 
combinations of local densities and currents up to second order in the derivatives. The method 
outlined in this Rapid Communication could be easily extended to functionals containing higher 
order terms, and could serve as a guide for constraining the time-odd part of the functional. 
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I. INTRODUCTION 



The density functional theory (DFT) {lj| has been very successfully employed in a wide 
variety of fields, from chemistry p] to condensed matter physics js-Sj. A somehow similar 
approach called the nuclear energy density functional (EDF) theory, has been developed in 
nuclear physics [6]. Even though the EDF theory is very often referred to as nuclear DFT, it 
is conceptually different due to the breaking of various symmetries in its current formulations 
(see e.g. [t-Q]). The EDF theory is a method of choice for describing the structure and the 
dynamics of medium-mass and heavy nuclei. Phenomenological nuclear functionals are now 
able to fit essentially all experimental nuclear mass data with rms deviations falling below 



0.6 MeV 



10|. Moreover the EDF theory is particularly well-suited for studying the dense 



inhomogeneous nuclear matter found in neutron stars and supernova cores. 

Nevertheless it has been known for a long time in condensed matter physics that the 
internal energy of a single electron as calculated by the DFT may be non-zero owing to a 
spurious interaction of the electron with itself (for a recent review, see for instance Section 
LB of Ref. lilj and references therein). This problem is not inherent to the DFT, but is 
induced by the use of approximate expressions for the electron exchange-correlation part 
which violate the Pauli exclusion principle. The presence of self-interaction causes various 
errors and is currently one of the main deficiencies of existing density functionals. For 
instance, this leads to an artificial spreading of the electron density owing to an incorrect 
asymptotic behavior of the Kohn-Sham potential. This spurious effect spoils the calculations 
of chemical properties. In particular, functionals contaminated by self-interactions fail to 
correctly describe the dissociation behavior of some homonuclear and heteronuclear diatomic 
radicals. The removal of self-interactions in the DFT has lead to the development of orbital- 



dependent density functionals 
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12|. 
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This self-interaction problem was also recognized a long time ago in nuclear physics 
14| . but it has been barely discussed in the literature until recently SQ-Q. Nuclear 



functionals have been traditionally derived from effective interactions using the Hartree-Fock 



approximation 
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20| which guarantees that in the one-particle limit, the internal energy 



correctly vanishes due to the exact cancellation between direct and exchange terms. On the 
other hand, the development of nuclear functionals which are not directly constructed from 
effective forces 2l|-|23l| calls for a closer examination of their internal consistency. 
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In this Rapid Communication, it will be shown how to eliminate self-interaction errors in 
the one-nucleon limit in semi-local nuclear energy density functionals containing all possible 
bilinear combinations of local densities, currents and their derivatives up to second order. 



II. SKYRME FORCE VERSUS NUCLEAR ENERGY DENSITY FUNCTIONAL 



Throughout this paper, pure nucleon states will be assumed. The more general formalism 



involving neutron-proton mixing has been developed in Ref . 2J| . The nuclear energy density 
functionals that we consider here are of the form 

E = -Ekin + -Ecoul + -Esky , (1) 

where i^kin is the kinetic energy, -Ecoui is the Coulomb energy and -Esky = / d'^r £^sky('') is 
the nuclear energy. The latter is a functional of the density matrix in coordinate space 
assumed to be purely diagonal in isospin space with elements Pq(r, a; r', a') with q = n,p for 
neutron and proton respectively (denoting the spin states by a, cr' = 1, — 1 for spin up and 
spin down respectively). In particular, Ss]iy{r) can be expressed in terms of the following 
local densities and currents (a is used to indicate the Pauli spin matrices and /i, u, n are 
indices of Cartesian tensors): 
(i) the density 

P^W = Pq{r,a;r,a) , (2) 



cr=±l 



(ii) the kinetic density 



r,ir)=J2 /ciV5(r-r')V-V'p,(r,a;r',a) , 



(3) 



(iii) the current-vector density 

(iv) the spin pseudovector density 

SqiT)= ^ Pq{r,a]r,o'){o'\a\a) , (5) 

(T,cr'=±l 
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(v) the spin pseudovector kinetic density 

TUr)= J2 f d'r'S{r-r')V-Vp,{r,a-y,a'){a'\a,\a) , (6) 

(vi) the spin current pseudotensor density 

J,,Ar)^-lYl [ d'r'S{r-r'){V,-V^)p,{r,ay,a'){a'\a^\a) , (7) 

a,a'=±l 

and (vii) the tensor-kinetic pseudovector density 

<t,<t'=±1 V 

It is convenient to introduce the isospin index i = 0, 1 for isoscalar and isovector quantities 
respectively. Isoscalar quantities (also written without any subscript) are sums over neutrons 
and protons (e.g. Po — P — Pn + Pp) while isovector quantities are differences between 
neutrons and protons (e.g. pi = pn — Pp)- The energy density £^sky constructed from all 
possible bilinear terms up to second order in the derivatives that are invariant under time 
reversal, space inversions and rotations, is given by 

^sky = 5^ (^r^^ + ) , (9a) 

t=0,l 



+ CfsfFt . (9c) 

The "time-even" part S^'^^^ ("time-odd" part £°^^) contains only even (odd) densities and 
currents with respect to time reversal. In situations for which the time-reversal symmetry 
is preserved, the time-odd densities jq{r), Sq{r), Tq{r) and Fq{r) must vanish. The couphng 

"constants" Cf and generally depend on the isoscalar density p = p„ + Pp as follows 

Cf = af + fcfp" , (10) 



Ct = al + hlp^ 



The spin-current vector density Jqir) appearing in Eq. ( l9b|) is defined by 



' qvK 



(12) 



where e^i/^ is the Levi-Civita tensor. Gauge (including Gahlean) invariance imposes the 



following relations 



26| 



(13) 



Historically the type of functionals given by Eqs. ([T]) and flUa|) -flUc |) was obtained from the 



Hartree-Fock approximation using effective zero-range interactions of the Skyrme type 
20, 
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27| 



1 1 

to(l + XQP,)5{r,j) + -ti(l + xiP,)^ [pI 5{rij) + 5{r,j) pi] 

+1(2(1 + X2Pa)^Pij ■ S{rij)pij + ^^3(1 + X3P^)p{r)°' 6(rij) 
+^Wo{ai+aj) -pij X 6{rij)pij 



+M [3(ai ■Pij)iaj -Pij) - {ai ■ aj)pl]S{rij) + 5{rij) [3{ai ■p^j){aj -Pij) - {ai ■ aj)pl] 



(14) 



where r^j = rj — r^, r = (rj + rj)/2, pij = — i/i(Vi — Vj)/2 is the relative momentum and 
Per is the two-body spin-exchange operator. The relations between the coupling constants 
in Eqs. fl9bl) and fl9cl) and the parameters of the effective force in Eq. f|T^ . can be found for 
instance in Table I of Ref. j24]. The parameters of the force are generally fitted to a selected 
set of nuclear data, mainly involving the time-even part of the functional. The coupling 
constants of the time-odd terms are then obtained from those of the time-even terms making 
use of the forementioned correlations (see e.g. Tables IV and VI in Ref. 26(]). However, since 
time-odd terms are not directly fitted to experimental data it is rather unlikely that nuclear 
properties involving those terms will be correctly described in this way. In particular, one 
of the main deficiencies of Skyrme forces is that they predict the occurence of spurious 
transitions to spin-ordered phases in the dense matter found in neutron stars and supernova 
cores 
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301 1 . It has been recently found that finite-size instabilities could arise at even 



lower densities thus spoiling self-consistent calculations in finite nuclei 
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32| 



Over the past decades, the conventional wisdom has been to regard the nuclear energy 
density functional as being more fundamental than effective forces and the development of 
microscopic functionals using many-body techniques with realistic nucleon-nucleon interac- 
tions has triggered a burst of activity |6|. Unfortunately such kinds of ab initio nuclear 
energy density functionals able to reproduce existing experimental nuclear data with the 
same degree of accuracy as effective interactions are not yet available. In the meantime 
semi-microscopic functionals have been constructed by first fitting the bulk part of the 
functional to realistic calculations of uniform infinite nuclear matter, and second adding a 
phenomenological surface part whose parameters are adjusted to reproduce properties of 



finite nuclei 
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33| . A complementary approach guided by effective field theories consists 



in constructing semi-local phenomenological functionals from combinations of local densities 
and currents up to a given order in the derivatives that are invariant under time reversal, 
space inversion, spatial rotations and gauge transformations 23 [. The expansion in terms 
of higher-order derivatives has been shown to converge rapidly 3^ . The coupling constants 
associated with each term are then directly fitted to a selected set of nuclear data {ssj. 
However, such kinds of unrestricted fits could lead to formal inconsistencies. In particular 
the resulting functional may allow a nucleon to interact with itself. As will be discussed in 
the next section, this unphysical feature can be avoided by imposing suitable constraints on 
the coupling constants of the functional. 



III. REMOVAL OF SELF-INTERACTION ERRORS 

The cancellation of spurious self-interactions in nuclear energy density functionals was 
studied in Ref. jl4|. The main assumption was that the nuclear functional can be calculated 
from the expectation value of some general Hamiltonian with a Slater determinant. However 
it remains to be proved that the exact functional can be obtained in this way. Moreover the 
nuclear functional was supposed to be purely local and depend on the densities ([2]) and the 
spin densities (|5]) only. In the following, the self-interaction problem will be reconsidered for 
any semi- local functionals given by Eqs. ([1]) and fl9a|) -fl9c |) . 

In the one-particle limit, the density matrix is simply given by pq{r,a;r' ,a') = 
where Lp^''\ra) is the single-nucleon wavefunction and the asterisk in- 
dicates complex conjugation. Inserting this expression into Eqs. (Elj-Q, substituting in 
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Eqs. ( pa|) - (pc|) and making use of the gauge invariance Eq. (fT3|) yield 



£7 even _i_ c 



odd 



Ap 



As\ 



1 2 



+ (CI - 2CI - ACi 



c 



(15) 



where cp denotes the two-component spinor wavefunction and the superscript f indicates 
the adjoint operation. Now the functional should be devoid of self-interactions whether 
the nucleon is free or is embedded in an external potential, i.e. for any single-nucleon 
wavefunction. The requirement that the nuclear energy -Esky vanishes thus leads to the four 
identities 



(16) 



Q + CI + C'q + Ci = A{C^' + Cf ^ + C^' + ci 



yAp 



tAs\ 



(17) 



Co" + CI - 2{C^ + Ci) - (C^ + Cf ) - A{Ct + Cf ^) = 



■iF\ 



(19) 



If the coefficients Cf and are allowed to depend on density according to Eqs. (ITOl) and 
(1 11 p . following the same analysis entails that Eq. flTBl) would have to be replaced by 



apQ + a{ + a'f^ + al = , 



(20a) 



K + K + bo + K = 



(20b) 



Note that these two conditions lead to Eq 



. (HM . but the reverse is not necessarily true. 



Equation flTB]) . was already noticed in Ref. [36|, even though its origin was not discussed. 

Equations f[TB]) - fll9l) can be automatically satisfied by first fitting the parameters of a 
Skyrme force given by Eq. (fT4|) and then calculating the corresponding coupling constants. 
But this usual fitting protocol imposes a lot more relations on the different terms of the 



functional than what is required by the cancellation of self-interactions (see e.g. Tables IV 
and VI in Ref. Actually these relations are generally not strictly enforced so that the 

functional may be contaminated by self-interactions. For example, the coupling constant C/ 



(hence also Cj due to gauge invariance) is set to zero in most Skyrme parametrizations 30 1. 
It has been suggested to also drop the isoscalar terms associated with Cq and C^^ in the 
time-odd part of the functional {20! . S^]. 

Self-interaction errors in the one-particle limit can contaminate systems consisting of 
many particles. Let us consider for instance cold fully polarized neutron matter. Assuming 
that all spins are aligned along the z-axis and neglecting the anisotropies induced by the po- 
larization, Eqs. (|^ -(|^ lead to the following expression for the energy density of polarized 
neutron matter 



ppol 
^NeuM 



n 



^ + {ci + CI + + + lic^ + ))p 

+{C^ + C'^+C^ + Ct)p' , ' (21) 



with the kinetic density 

rPor=hQ^2Y/3^,/3 (22) 



It can thus be seen from Eqs. (|T6|) and fl2T|) that the error caused by self-interactions is given 
by 

S£^tM = iC'o+C',+C^ + Ct)p' . (23) 

As expected the magnitude of this spurious term grows with increasing density and could 
thus lead to large deviations when such functional are applied to the high-density matter 
found in neutron stars. In particular if ^i^^eliM < 0' self-interaction errors will drive a 
ferromagnetic collapse of neutron matter with the energy fl2T]) decreasing asymptotically 
with increasing density. Such a pathological behavior is obviously contradicted by neutron- 
star observations. 



IV. CONCLUSIONS 



Phenomenological nuclear functionals based on Skyrme forces have been very successful 
in describing various nuclear systems, from finite nuclei to neutron stars and supernova 
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cores. Nevertheless the use of effective forces introduces tight correlations between different 
terms of the functional which hamper further improvements. For instance, adding density- 
and momentum-dependent terms in ti and ^2 contributes tojthe eight coupling constants 
Cl, Cj, C^^ and C^^ , and also generates many new terms [28]. Improving one part of the 
functional may therefore deteriorate other parts. 

A different approach to improve Skyrme functionals systematically consists in construct- 
ing semi-local functionals from all possible combinations of local densities and currents up 
to a given order in the derivatives, that are invariant under time reversal, space inversion, 
spatial rotations and gauge transformations [23]. The coupling constants appearing in the 
functional are treated as free adjustable parameters and are directlty fitted to a selected 
set of nuclear data. Unfortunately such functionals inspired by effective field theories may 
allow a single nucleon to interact with itself thus indicating a violation of the Pauli exclusion 
principle. This unphysical feature can be avoided by imposing suitable restrictions on the 
coupling constants. These constraints have been derived for semi-local functionals up to 
second order in the derivatives and are given by Eqs. ( fT6l) -( |T9l) . The method outlined in 
this Rapid Communication can be easily extended to functionals of higher order. 

Self-interaction errors in the one-particle limit can also contaminate systems consisting 
of many particles. For instance, it has been shown that in some cases these errors can lead 
to a catastrophic ferromagnetic collapse of neutron matter in contradiction with neutron- 
star observations. More generally, because a single nucleon breaks time-reversal symmetry, 
preventing a nucleon from interacting with itself necessarily requires the presence of time- 
odd densities and currents in the functional. The cancellation of self-interactions could thus 
serve as a guide for constraining the time-odd part of the functional, which still remains 
poorly known. 

Even though Eqs. (fT6l) - (fT9l) ensure the vanishing of the internal energy in the limit of 
a single nucleon, other kinds of self-interactions and self-pairing could still arise in many- 
body calculations Nuclear energy density functionals obtained from the Hartree-Fock 
method with effective forces may also be spoiled by such many-body self-interactions. In 
particular, it is well-known that the density-dependent term in the Skyrme force f|T^ is 
formally inconsistent from the many-body theory point of view. Whereas its contribution to 
the ground-state energy vanishes in the one-nucleon limit, this term contains spurious many- 
body self-interactions which can become manifest when dealing with more than one particle. 
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Correcting for these many-body self-interactions would provide additional constraints on 
nuclear energy density functionals. 

Acknowledgments. This work was financially supported by FNRS (Belgium), Commu- 
naute frangaise de Belgique (Actions de Recherche Concertees) and by CompStar (a Re- 
search Networking Programme of the European Science Foundation) . The author is grateful 
to J.M. Pearson for valuable comments. 



[1] W. Kohn, Rev. Mod. Phys. 71, 1253 (1999). 

[2] W. Koch and M.C. Holthausen, A Chemist's Guide to Density Functional Theory, 2nd ed. 

(Wiley- VCH, Heidelberg, 2001). 
[3] R. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University 

Press, Cambridge, 2004). 
[4] A. Bulgac, Phys. Rev. A 76, 040502(R) (2007). 
[5] G. Rupak and T. Schafer, Nucl. Phys. A 816, 52 (2009). 

[6] J. E. Drut, R. J. Furnstahl, L. Platter, Prog. Part. Nucl. Phys. 64, 120 (2010). 
[7] J. Engel, Phys. Rev. C 75, 014306 (2007). 

[8] J. Messud, M. Bender, E. Suraud, Phys. Rev. C 80, 054314 (2009). 

[9] T. Duguet and J. Sadoudi, J. Phys. G 37, 064009 (2010). 
[10] S. Goriely N. Chamel, and J. M. Pearson, Phys. Rev. C 82, 035804 (2010). 
[11] S. Kiimmel and L. Kronik, Rev. Mod. Phys. 80, 3 (2008). 
[12] J.P. Perdew, A. Zunger, Phys. Rev. B 23, 5048 (1981). 

[13] G. Ripka, Lectures on nuclear models with density dependent effective interactions, preprint 
DPh-T/69-54. In: Service de Physique Theorique, Centre d'Etudes Nucleaires de Saclay 
(1969). 

[14] S. Stringari and D.M. Brink, Nucl. Phys. A304, 307 (1978). 
[15] M. Bender, P.-H. Heenen and P.-G. Reinhard, Rev. Mod. Phys.75, 121 (2003). 
[16] D. Lacroix, T. Duguet and M. Bender, Phys. Rev. C 79, 044318 (2009). 
[17] M. Bender, T. Duguet and D. Lacroix, Phys. Rev. C 79, 044319 (2009). 
[18] T. Duguet, M. Bender, K. Bennaceur, D. Lacroix and T. Lesinski, Phys. Rev. C 79, 044320 
(2009). 



10 



[19] J. Erler, P. Kliipfel and P.G. Reinhard, J. Phys. G 37, 064001 (2010). 
[20] J. R. Stone and P.G. Reinhard, Prog. Part. Nucl. Phys. 58, 587 (2007). 
[21] S.A. Fayans, S.V. Tolokonnikov, E.L. Trykov, D. Zawischa, Nucl. Phys. A676 (2000) 49. 
[22] M. Baldo, L. Robledo, P. Schuck and X. Vinas, J. Phys. G37, 064015 (2010). 
[23] B.G. Carlsson, J. Dobaczewski, M. Kortelainen, Phys. Rev. C78, 044326 (2008). 
[24] E. PerUnska, S.G. Rohozinski, J. Dobaczewski, W. Nazarewicz, Phys. Rev. C69, 014316 
(2004). 

[25] Y. M. Engel, D. M. Brink, K. Goeke, S. J. Krieger, and D. Vautherin, Nucl. Phys. A249, 215 
(1975). 

[26] J. Dobaczewski and J. Dudek, Phys. Rev. C 52, 1827 (1995); 55, 3177(E) (1997). 
[27] T. Lesinski, M. Bender, K. Bennaceur, T. Duguet and J. Meyer, Phys. Rev. C 76, 014312 
(2007). 

[28] N. Chamel, S. Goriely, and J. M. Pearson, Phys. Rev. C 80, 065804 (2009). 
[29] Li-Gang Cao, G. Colo, H. Sagawa, Phys. Rev. C 81, 044302 (2010). 
[30] N. Chamel and S. Goriely, Phys. Rev. C 82, 045804 (2010). 

[31] T. Lesinski, K. Bennaceur, T. Duguet, J. Meyer, Phys. Rev. C 74, 044315 (2006). 
[32] M. Kortelainen and T. Lesinski, J. Phys. G 37, 064039 (2010). 

[33] L. G. Cao, U. Lombardo, C. W. Shcn and N.V. Giai, Phys. Rev. C 73, 014313 (2006). 

[34] B.G. Carlsson and J. Dobaczewski, Phys. Rev. Lett. 105, 122501 (2010). 

[35] M. Kortelainen, T. Lesinski, J. More, W. Nazarewicz, J. Sarich, N. Schunck, M. V. Stoitsov, 

S. Wild, Phys. Rev. C82, 024313 (2010). 
[36] M. Bender, J. Dobaczewski, J. Engel, and W. Nazarewicz, Phys. Rev. C 65, 054322 (2002). 
[37] P. Kliipfel, P.-G. Reinhard, T. J. Biirvenich and J. A. Maruhn, Phys. Rev. C 79, 034310 

(2009). 



11 



